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We generalise a non-classicality test described by Kot et al. [Phys. Rev. Lett. 108, 233601 
(2010)], which can be used to rule out any classical description of a physical system. The test is based 
on measurements of quadrature operators and works by proving a contradiction with the classical 
description in terms of a probability distribution in phase space. As opposed to the previous work, 
we generalise the test to include states without rotational symmetry in phase space. Furthermore, 
we compare the performance of the non-classicality test with classical tomography methods based 
on the inverse Radon transform, which can also be used to establish the quantum nature of a 
physical system. In particular, we consider a non-classicality test based on the so-called filtered 
back-projection formula. We show that the general non-classicality test is conceptually simpler, 
requires less assumptions on the system and is statistically more reliable than the tests based on the 
filtered back-projection formula. As a specific example, we derive the optimal test for a quadrature 
squeezed single photon state and show that the efficiency of the test does not change with the degree 
of squeezing. 
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I. INTRODUCTION 

Nearly a century after its foundation [1, 2], Quantum 
Mechanics is a well established theory for microscopic and 
mesoscopic phenomena. Nonetheless, the logical struc¬ 
ture of quantum mechanics is still challenging our per¬ 
ception of the theory, spurring investigations of a deeper 
understanding of the quantum theory and its relation to 
the classical perception of reality [3, 4]. In particular, the 
conceptual differences between the quantum and classical 
world remain a fundamental topic of investigation, in part 
due to the potential applications of quantum phenomena 
for quantum information technology [5]. To this end, ex¬ 
periments have sought to demonstrate effects which are 
truly non-classical , but the challenge in this case is how 
to prove that something really is non-classical. Clearly, 
the mere consistency with predictions from quantum the¬ 
ory is not sufficient to claim the observation of genuine 
quantum effects as this does not rule out a classical de¬ 
scription of the same phenomena, although the classical 
explanation may seem different from the quantum me¬ 
chanical. On the contrary, throughout this article, we 
shall use the definition that an experimental observation 
of a genuine non-classical effect should only assume con¬ 
cepts from classical physics and then demonstrate that 
these assumptions lead to a contradiction with the ob¬ 
served experimental result. In other words, the goal of 
the non-classicality test is to convince a physicist trained 
in classical physics, but completely ignorant of quantum 
theory, that his/her perception of the world is incorrect. 
A non-classicality test fulfilling this requirement based 
on measurements of quadratures operators, e.g., position 
x and momentum p , was presented in Ref. [6] . Further¬ 


more, the test was used to demonstrate that the electric 
and magnetic field of a single photon state, cannot be 
described classically. In this article, we generalise the 
approach of Ref. [6] and develop the theory for genuine 
non-classicality tests based on measurements of quadra¬ 
tures operators for a physical system. As opposed to 
Ref. [6], which only considered rotationally symmetric 
states, we consider states of arbitrary shape. Further¬ 
more, we compare the obtained results with more con¬ 
ventional approaches relying on tomographic reconstruc¬ 
tions, e.g. inverse Radon transformations, and demon¬ 
strate that the tests developed here are both conceptu¬ 
ally simpler, require fewer assumptions on the system, 
and are statistically more reliable. 

One of the most impressive examples of a genuine non- 
classical effect is the violation of Bell’s inequality. In 
Bell’s inequality, two simple classical assumptions are 
made, locality and realism, and this is shown to lead to 
a contradiction with the observed experimental results. 
Remarkably, this can be done without any assumptions 
about the physics of the underlying physical system. On 
the other hand, Bell’s inequality is often quite difficult 
to violate in practice [7, 8]. We, therefore, take a differ¬ 
ent approach and consider non-classicality tests based on 
measurements of quadrature quantities of the form 

Qg = cos 0 a; + sin0p (1) 

for a single system. We shall assume that Qg exist for 
any 6 £ [0, 7r] , in the sense that they are well defined 
measurable observables of the physical system. In par¬ 
ticular, we are assuming that there exist two independent 
observables x and p, which might be canonical variables 
although this is not necessary. As considered in Ref. [6], 
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this is for example the situation for optical fields mea¬ 
sured by homodyne detection. In this case, a completely 
classical analysis reveals that Qg are the measured ob¬ 
servables with the angle 9 controlled by the phase of an 
interfering laser field and x and p are, e.g., the electric 
and magnetic field components of a single mode of the 
measured optical field. Another example can be found in 
the emerging field of opto-mechanics [9, 10] where a lot 
of effort is devoted to proving that mechanical oscillators 
can have non-classical behaviour. Here, a measurement 
of the position x m (t) at a time t will have the form in 
Eq. (1) with 9 = co m t where Lo m is the vibration frequency 
of the oscillator and x and p are the (suitably rescaled) 
position and momentum operators at time t = 0. Finally, 
operators of the form in Eq. (1) can also be measured for 
atomic ensembles where x and p represents two differ¬ 
ent spin components of the ensemble [11]. In principle, 
operators of the form in Eq. (1) can be used to show a 
violation of Bell’s inequality [12] but this requires prepar¬ 
ing quantum states, which are highly difficult to produce. 
To avoid this, we make a stronger assumption compared 
to Bell’s inequalities and assume that the measured ob¬ 
servables are described by the relation in Eq. (1) as pre¬ 
scribed by the underlying physics of the system. This 
assumption makes our non-classicality inherently weaker 
than Bell’s inequality, but on the other hand enables the 
violation of the non-classicality test under less stringent 
experimental requirements. 

As a second assumption underlying any classical de¬ 
scription of the system, we assume that the quadratures 
defined by Eq. (1) commute with one another. As a 
result, these two hypotheses give the complete picture 
of the observable quantities that we will need for our 
non-classicality test. Mathematically speaking, it means 
that we postulate the algebra of observables to be the 
commutative algebra generated by the rotated quadra¬ 
tures. These hypotheses underlie the classical descrip¬ 
tion of the system, which we will be able to rule out 
by deriving a contradiction from experimental data. In 
this description where x and p are commuting quantities, 
the state of the system will be characterized by a joint 
^-probability distribution, also called the phase space 
distribution W c (x,p). Then, the contradiction consists 
in showing that such a distribution does not exist for the 
measured experimental outcomes. 

In quantum mechanics, it is well known that there does 
not exist a joint probability distribution for x and p since 
the position and momentum of a particle do not com¬ 
mute. The closest thing that one can have is the Wigner 
function W(x,p) [13], which correctly predicts the proba¬ 
bility outcome of measurements of the observables Qg but 
which cannot be considered a proper distribution func¬ 
tion since it can attain negative values (see Fig. 1). The 
negativity of the Wigner function will be the underlying 
quantum mechanical reason for the failure of the classi¬ 
cal description in the non-classicality test derived below. 
In quantum mechanics, the Wigner function is, however, 
not the only possible phase space distribution and nu¬ 


merous others have been considered as the basis for non- 
classicality criteria. For instance the P-function repre¬ 
sents the expansion of a field on coherent states and is 
thus negative for any state with an uncertainty squeezed 
below the level of a coherent state. Since, according to 
quantum theory, classical light sources produce coherent 
states, any state produced by classical sources will have a 
positive P-function. As a consequence, states with neg¬ 
ative P-functions are often referred to as non-classical 
and based on this definition numerous observable non- 
classicality criteria have been derived [14-17]. These cri¬ 
teria, however, inherently rely on quantum mechanics 
since the underlying criteria is based on having smaller 
uncertainty than expected from the quantum mechani¬ 
cal uncertainty relation. With the definition used here, 
i.e., that a non-classicality criteria should only rely on 
concepts from classical physics, tests based on the nega¬ 
tivity of the P-function are thus not applicable. From a 
mathematical perspective, however, the test that we will 
derive is highly related. In particular our test is essen¬ 
tially the same as the one derived for the P-function in 
Refs. [ 1, 17], but it is translated to Wigner functions in 
such a way that quantum mechanics is not assumed, thus 
leading to a stronger test. 

II. SIMPLE NON-CLASSICALITY TEST 
A. Premise 

Before proceeding, we shortly review the non- 
classicality test presented in Ref. [6], which is based on 
previous work by Bednoraz and Belzig [18]. We will em¬ 
phasise the fundamental ideas behind the test and outline 
how it can be generalized into the more general test in¬ 
troduced below. 

As already discussed in the introduction, a genuine 
non-classicality criterion should directly disprove the 
classical description without logically relying on quantum 
mechanics. Because the state of the system is classically 
characterized by a probability distribution W c {x,p) over 
the phase space of x and p, the expectation value of any 
function gf(x,p) is 

{&/) = J j dxdpW c {x,p)srf(x,p). (2) 

Mathematically speaking, Eq. (2) is the integral repre¬ 
sentation of positive, normalized, linear functionals on 
the <7*—algebra generated by the observables. Accord¬ 
ingly, any non-negative quantity JF(x,p) > 0 should have 
a non-negative expectation value, 

{&) = J J dxdpW c {x,p)&(x,p) > 0, (3) 

because the integrand is non-negative, W c (x,p) being a 
probability distribution. The violation of Eq. (3) dis¬ 
proves W c (x, p) as a probability distribution and thus 
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FIG. 1 (colour online), (a) Wigner function W(x,p) of the single photon state of a single mode of light. It is not a probability 
distribution as it should be in the classical framework, since it takes on negative values, (b) Top view of the same function. 
The colour gradient allows to identify the negative region in the centre. The cut represented by the line is the quadrature 
Qg. The marginal distribution corresponding to the probability distribution for measurement outcomes for Qg is obtained by 
integration of W(x,p) along the perpendicular directions to the cut, as shown by the arrows. 


rules out any classical description of the system. We 
cannot, however, directly determine (&) since we cannot 
measure x and p simultaneously and we have to infer it 
from measurements of Qg. 

In Ref. [6], the test function & was chosen to be the 
square of a polynomial in the squared radius r 2 = x 2 +p 2 
since the underlying physical state was assumed to be 
rotationally symmetric. When the polynomial in r was 
expanded, the various moment of r could be expressed in 
moments of the rotated quadratures Qg. Consequently, 
the mean value (J^) together with its experimental er¬ 
ror were directly computable from the experimental out¬ 
comes. In Ref. [6], the developed test was applied to 
the experimental data for an approximate single photon 
state, attaining the violation of Ecp (3) by nineteen stan¬ 
dard deviations [6]. 

B. General formulation 

We will now generalise the test function considered in 
Ref. [6] so that it does not rely on rotational symmetry. 
Consequently, the generalized formulation can exploit the 
freedom in the choice of phases for the quadratures to 
construct an optimal test function for an arbitrary quan¬ 
tum state. 

Let us define the problem in more rigorous terms. To 
formulate a test, we need a suitable class of non-negative 
test functions & defined on the phase space and a method 
to determine (J^), together with its empirical error, di¬ 
rectly from the measured quadrature moments. The gen¬ 
eral class of smooth non-negative functions on the phase 
space is too large for our purposes and we therefore re¬ 
strict the test to non-negative polynomials in the vari¬ 


ables x and p. Note, that such polynomials can be peaked 
in a region around the negativity of WQc^p) and since 
W {x, p) typically vanishes faster than any power law out¬ 
side its significant support, the divergent terms of the 
polynomials are irrelevant (see Fig. 2). Consequently, we 
do not loose much generality by restricting to polyno¬ 
mials, and as we will show below, it allows us to make 
an explicit representation of the test function in terms of 
quadrature quantities. 

Motivated by these considerations, we consider test 
functions of the form 

N i \ 2 

1 + ^^VV ■ (4) 

i=1 k =0 / 

Here we have, for simplicity, put the zeroth order coeffi¬ 
cient equal to one, since it is irrelevant to the effectiveness 
of the test as discussed below. Furthermore, we choose to 
ensure the positivity by only considering squared poly¬ 
nomials. This choice simplifies the test as discussed in 
Appendix A, but does not exhaust all the possible, non¬ 
negative polynomials. This is related to the 17 th Hilbert 
problem [17, 19]. However, since a polynomial of suffi¬ 
ciently high order can approximate an analytical function 
sharply peaked at any particular value where W(x,p) is 
negative, we believe that considering this restricted set 
does not degrade the efficiency of the test significantly. 

To have an experimentally applicable test, J^( x,p) 
needs to be expressed in terms of powers of measurable 
quantities. For a fixed i, Q' g directly gives a sum of (i + 1) 
moments of the form x l ~ k p k with k = 0, ..., i. Hence, if 
one measures m > i+1 different angles 6j, one can always 
combine those linear equations by simple linear algebra 
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to arrive at 


and using Eq. (5), we get 


m 

x i- kp k = J2 T i . kJ Q i g ., (5) 

j=i 

where T^kj are real-valued coefficients which have to sat¬ 
isfy the linear constraint 

^2 ( k ) cosI tanfc % = 4,fc' Vi. (6) 

Note that at least (i+1) linearly independent quadratures 
are needed to express all (i + 1) moments of x and p, 
but if m > i + 1, the constraint does not give a unique 
solution and there is a freedom in the choice of the Ti-k.j 
coefficients. 

Before proceeding, we stress that we have assumed 
commutativity of x and p in attaining Eq. (5), i.e., we 
assume a classical description of the quadratures in terms 
of x and p. From Eq. (5), we can express the test func¬ 
tion in terms of the quadratures Q l s , which are directly 
measurable by assumption. Accordingly, we have explic¬ 
itly shown that the test function ^{x,p) is in the algebra 
of observables and through Eq. (5), we are able to com¬ 
pute its mean value. Finally, it is worth noting that this 
strikingly simple linear relation is the very heart of the 
procedure we consider. This is the trick that allows us to 
avoid the complicated machinery of the inverse Radon al¬ 
gorithms for the inverse problem of finding W(x,p) from 
the measured quadratures Qg as we will discuss in Sec. 
III. 

In summary, once the polynomial order of the test func¬ 
tion is fixed to 2TV, the test is characterized by three sets 
of parameters: 

• the set of free coefficients that define the test 
function over the phase space, 

• the set of m > 27V + 1 free phases {6j} which iden¬ 
tify the quadratures we measure, 

• the set of coefficients T^j , which allow to express 
the x l ~ k p k moments in terms of Qg. . These are not 
completely free to vary but has to satisfy the linear 
constraint in Eq. (6). 

We will formally refer to a test, T, as the string of pa¬ 
rameters 


T = (TV, Di-k , 9j, Ti-kj). 


( 7 ) 


2N i 

^(x, P ) = i+Y,J2 c ^ xl ~ k p k 

i— 1 k =0 

m 2N / i 

= 1 +EEqUE^ t ^ 

j =1 i =1 \k—0 

m 

= 1 + E^(Q* J ). (8) 

i=i 

For simplicity, we have here introduced the coefficients 
fV ?;:k = 2 Di-k T — l - 7b ;:h Hq ’ : ly ', where Di-k 

vanishes for i > N. Furthermore, we have introduced the 
set of m polynomials, [Qg ), (one for each phase cut) 
defined as 


i =1 

Eq. (8) shows that according to the relation in Eq. (5), 
the test function can be expressed as a polynomial of 
quadratures, i.e., as a sum of 27V + 1 polynomials of de¬ 
gree 27V, one for each independent phase cut. 

To evaluate the expectation value of the final expres¬ 
sion in Eq. (8), we consider each quadrature as an in¬ 
dependent random variable distributed according to its 
marginal. In principle, if one has a very large number of 
cuts (in particular if the number of cuts m exceeds the 
minimal required value m > 27V + 1), one could argue 
that nearby cuts will have similar moments and this may 
be used in a statistical analysis to reduce the variance on 
the experimental data. On the other hand, we prefer to 
make as few assumptions as possible and therefore con¬ 
sider each cut to be an independent variable, distributed 
according to its marginal. Thus, for fixed phase cut, the 
expectation value Qg)) is formally given by integra¬ 
tion with the marginal distribution and is experimentally 
estimated as 

wm) * Enil (10) 

with \Qg\ being the n-th outcome within the dataset 
consisting of a total of Mj measurements of the random 
variable Qg. Accordingly, we can compute the average 
value of the test function & as 

m 

(JT) = 1 + £(■*? (<&,)>, (ii) 

3 =1 


2 N 

■*3(^) = E<4 



( 9 ) 


We will now describe, in more details, how to express 
the expectation value of the test function and its relative 
empirical error in terms of the measurable quadratures 
Qg .. Expanding the square in the test function in Eq. (4) 


which according to a classical description, must be posi¬ 
tive (See Eq. (3)). 

Quantum mechanically, (&) is given by Eq. (2) with 
W c (x,p) being the Wigner function W{x,p) describing 
the phase-space distribution. For a general density ma¬ 
trix, p, the Wigner function is defined as 

1 r°° 

W(x,p) = —J (x + y\p\x-y)e~ 2lpv,h dy, (12) 
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where \x ± y) are eigenstates of the position operator 
with eigenvalue x ± y. The Wigner function is often re¬ 
ferred to as a quasi-probability distribution since it can 
attain negative values and is thus not a true probability 
distribution. This possibility of the Wigner function be¬ 
ing negative is the key aspect of our non-classicality test. 
With a negative Wigner function, Eq. (2) and the equiv¬ 
alent form in Eq. (11) are no longer guaranteed to be 
positive leading to contradiction with classical physics. 
Given the Wigner function describing a quantum sys¬ 
tem, we can find the set of coefficients that mini¬ 
mizes (jF), by solving a linear optimization problem as 
described in Ref. [6]. In an experimental test, however, 
we also have to consider the empirical error of our esti¬ 
mate of (^). Since each random variable, Qg , is treated 
as statistically independent for different 0, we can esti¬ 
mate the variance of (J?) as 


Note, however, that even though the optimization as¬ 
sumes the quantum mechanical description of the state, 
the test itself is still completely classical, and the assump¬ 
tion of quantum mechanics at this stage does not decrease 
the strength of the test. In an experiment, the optimal 
choice of phase cuts has to be determined before the ex¬ 
perimental measurements while the optimization of the 
other free parameters can be performed afterwards and 
may be optimized over the experimental data. If this 
is done, the statistical analysis of the experiments should 
also include the freedom in varying these parameters. For 
simplicity, we will just assume that one uses the coeffi¬ 
cients expected from theory, in which case we do not have 
to worry about this complication. 

C. Rotationally invariant states 


m 1 

(«(<?*)> - mQe)) 2 ), ( 13 ) 

3 =1 C 


where we have assumed for simplicity that the same num¬ 
ber of measurements M c is performed for each cut. This 
is sensible because of the statistical independence as¬ 
sumed, which means that no bias is justified (for a formal 
justification of this choice see Appendix A). 

We are interested in classifying how well the various 
non-classicality tests perform. To do so, we introduce 
the following quantity for any non-classicality test, 


S[T] 


V (&) 

lim - , 

M ->°o o&y/M — TO 


(14) 


which depends on all the characteristics of the particu¬ 
lar test T, but not on the total number of measurements 
M = to. M c because a & scales asymptotically as 1/ y[M. 
In fact, as shown in Appendix A, a jr goes exactly as 
1 /y/M — m for the optimal test and the M dependence 
thus disappears even before taking the limit. Note that 
because Q involves a ratio, any overall proportionality 
constant in the test function is irrelevant. We have there¬ 
fore chosen to set the first coefficient in the polynomial 
in Eq. (4) equal to unity to get rid of any overall scaling 
factor. 

The optimal test among all possible tests T (See 
Eq. (7)) is defined as the one, which minimizes Q in 
Eq. (14). Formally, the search for the optimal test in¬ 
volves an algorithm to choose a suitable test among a 
class of logically equivalent tests. Such an optimization 
is generally hard to do analytically and we will there¬ 
fore treat it numerically. Since we want to prove the 
non-classicality of the system, it is natural to optimize 
the test using the quantum mechanical description of its 
state, i.e., by assuming the suitable Wigner function to 
be the correct description of the unknown distribution 
W(x,p). We can then employ simple symmetry consid¬ 
erations about W(x,p) in the optimization, as we show 
below for rotational invariant states and squeezed states. 


If the distribution W(x,p), describing the state, is ro¬ 
tationally invariant, we naively expect that all the pa¬ 
rameters of the optimal test function should exhibit this 
rotationally invariance, e.g., that the optimal phase cuts 
{0j} should be uniformly distributed. This was assumed 
in Ref. [6] but not formally shown. We will now show 
that this is indeed the case. 

We explicitly define the rotational version of the test 
function as 




W 2J 


i+ £ d i' 


2 i 


(15) 


with r 2 = x 2 + p 2 being the squared radius in the phase 
plane. Apart from the fact that we have formally reduced 
the set of coefficients D r k to the lighter set of di, we make 
a further reduction when we relate the radial moments 
to the quadrature moments in a similar way as done in 
Eq. (5). We seek a linear relation of the form 

m 

7-21 = £*i;jQe’> ( 16 ) 

3=0 

with the set of coefficients, t rj , being the analogue of the 
set of Ti-kj in Eq. (5). The existence of such a linear 
relation, provided a sufficient number of phase cuts, can 
be deduced from arguments similar to the ones leading 
to the relation in Eq. (5). Note that odd powers of the 
quadratures do not appear in the above formula, since the 
parity symmetry of r 2 * guarantees the null contribution 
of such terms. Denoting rotationally invariant tests as 
Trad? we prove in Appendix A the following theorem. 

Theorem 1. Whenever W(x,p) is rotationally symmet¬ 
ric, and takes on negative values, the optimal test, i.e., 
the test T that minimizes Q, is a radial test T ra d and the 
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cuts are chosen in such a way that 

dj = ^ 3 = 0, -,m- 1 (17) 

hi = u= ( 2 - ) ^ Vj. (18) 

Furthermore Q is independent of the number of cuts pro¬ 
vided that m> N + 1 . 

The theorem states that the optimal test function is 
a radial polynomial and the cuts should be distributed 
uniformly over the phase space whenever W(x,p) is ro- 
tationally symmetric. Note that the numerical optimiza¬ 
tion of the test is a lot simpler than in the general case 
because the optimal test function can be put into the 
form in Eq. (15), which has fewer free parameters. Yet, 
the reader should bear in mind that with the choice given 
by Eqs. (17) and (18), the linear relation in Eq. (16) 
holds and can be applied even for non-symmetric states. 
Hence, the assumptions made can never lead to a state 
incorrectly being assigned as non-classical even if the as¬ 
sumptions are not fulfilled. Applying a rotationally in¬ 
variant test to a non rotationally invariant state will thus 
only lead to a less efficient test but not to an incorrect 
result. 

For the test T ro d, we can express the test function, 
using Eqs. (16) and (18), as 

m 

^ = l + J2^(Qo J ), (19) 

i=o 

where Jif is now the same polynomial for all the phases, 
contrary to the general case in Eq. (11). Defining c, = 
2 di + X] a +a'=i d a d a 'i we can explicitly express this poly¬ 
nomial as 

2N 

Jt?{Q e )=Y, c ^Q 2 e- ( 2 °) 

z — 1 

Otherwise all considerations about the evaluation of Q 
are the same as for the general test. Most importantly, 
we treat every cut as statistically independent, even in 
the case of actual rotational symmetry, because we do 
not want to assume this symmetry. 

For the single photon state considered in Ref. [6], which 
corresponds to the first excited state of a harmonic oscil¬ 
lator, the Wigner function is 

W(x,p) = ^ (%(x 2 +p 2 ) - l) exp [-(a; 2 +p 2 )\■ (21) 

We have numerically optimized Q with respect to T for in¬ 
creasing polynomial order. Fig. 2 shows the optimal Q as 
a function of N, the order of the test function being 2N. 
We find a violation of Eq. (3) for N > 4 in agreement with 
what was found in Ref. [6]. Furthermore, we see that the 
graph plotted in Fig. 2 exhibits an asymptotic saturation 
of Q. Accordingly, in the case of arbitrarily high number 



FIG. 2 (colour online). Violation of non-classicality for the 
first excited state of a harmonic oscillator, e.g, a single pho¬ 
ton state. The figure show the optimal expectation value 
of the test function relative to its standard deviation ( Q ) 
as a function of the order of the polynomial in Eq. (15). 
Experimentally this may result in a violation of classicality 
by — vM — mQ standard deviations, with M being the to¬ 
tal number of measurements and m > N + 1 the number of 
measured quadratures. 

of measurements, increasing the polynomial degree be¬ 
yond N ~ 16 gives more degrees of freedom for the test 
but does not decrease Q significantly. This demonstrates 
that polynomial test functions are efficient. In principle, 
more general, bounded and smooth test functions could 
be employed but the results obtained here suggest that 
this is not necessary. 

While Q quantifies the relation between different tests, 
it does not directly give the statistical strength, (&) /o& 
of a test. The latter can, however, be obtained by sim¬ 
ply multiplying Q with \JM — m. As stated above Q is 
independent of the number of cuts, but because of the 
dependence of m in \JM — in it will in general be bet¬ 
ter to have the least number of cuts, be., m = N + 1, 
although this dependency will be small if a large num¬ 
ber of measurement are performed M m. In order 
to keep m as small as possible we thus find from Fig. 2 
that is minimised by working around N = 16 

although higher order tests IV > 16 will produce essen¬ 
tially similar results if sufficiently many measurements 
are performed M m. For imperfectly prepared states, 
however, the area in phase space where the Wigner func¬ 
tion is negative is smaller. Hence in this case it may be 
desirable to work at a higher N where the function can 
be tailored to be more narrow [6]. In Fig. 3, we plot the 
optimal test function for the minimum argument N = 16. 
It is seen that the test function is peaked at the negative 
region of the Wigner function and that the tails of the 
polynomial are outside the significant support of W(x,p) 
and are thus not affecting the test. 
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FIG. 3 (colour online). Radial profile of W(x,p) for the single 
photon state (red dashed curve) and of the optimal test func¬ 
tion IF{x,p ) (light blue solid curve) for N = 16. The inset 
shows the corresponding 3D plot of x,p ). 


D. Squeezed states 

We now turn to the case where W{x,p) is not rota- 
tionally symmetric. As a specific example, we consider a 
squeezed single photon state and show that the optimal 
test function differs significantly from a rotationally in¬ 
variant test. Let W(x,p) be the Wigner function of the 
single photon state. The Wigner function, W\(x,p) of 
the squeezed state is given by the rescaling 

W\ (x, p) = W{ Ax, p/X ), (22) 

where A characterizes the amount of squeezing in either 
x (A > 1) or p (0 < A < 1). The unsqueezed state corre¬ 
sponds to A = 1. Fig. 4 shows the result of the squeezing 
transformation; the Wigner function of the single photon 
state is stretched into an ellipse for the squeezed state. 

The squeezed state still exhibit elliptical symmetry but 
the marginal distributions are no longer identical as they 
were for the rotationally invariant state. The extra free¬ 
dom in the choice of cuts present in the general test is 
therefore of importance. As an ansatz for a test function, 
we consider 


FIG. 4 (colour online). Sketch (not to scale) of W(x,p) (left) 
and W\(x,p) (right), together with a representation of the 
optimal distribution of phase cuts. The cuts are uniformly 
distributed for the single photon state and stretched along 
the x axis for the squeezed state. Note that the cuts stretch 
along the opposite direction of stretching for W\(x,p). This 
reflects that most cuts are placed at angles where the distri¬ 
bution varies the fastest. The inner/blue region is where the 
distribution takes on negative values. 


p —I p/X such that 


rjf = ^ ti-j (cos 9j Xx + sin 6j p/X ) 

3=0 


2 i 


= ^ ti-j (Acos 9j) 

3=0 


2 i { tan 0j 

x H- 


A 2 


2 i 




(24) 


where 9 £ [0,7r) satisfies tan (A, =* taaOj . geometrically, 
this transformation 9 —> 9 is the polar representation of 
the rescaling that we have performed on the phase plane. 

In Appendix A, we show that not only is the average 
value of independent of A, but also the same is true 
for the figure of merit Q. This is formally stated in the 
following theorem. 


&\{x,p) 


[N/2] 


1+ d, 


\ \ ) 


(23) 


with r 2 = (Aa;) 2 + {p/X) 2 being the rescaled radius in the 
phase plane. Eq. (23) is a rescaling of the radial functions 
in Eq. (15) by the squeezing parameter, A . Note that 
with this type of test function, the average value of 
with a phase space distribution W\ is independent of A. 
The polynomial expression of &\(x,p) with respect to 
the rotated quadratures is obtained in a straightforward 
manner from Eq. (16) with the substitution x —> Xx and 


Theorem 2. (Invariance under squeezing) The mini¬ 
mum of Q for the squeezed state is independent of the 
squeezing parameter A and equal to the minimum for the 
rotationally symmetric state. The optimal test function 
is given by Eq. (23) with di being the same as the optimal 
ones for the single photon state. Furthermore 9j and t^j 
given in Eqs. (17,18) transform to the “rescaled’' values 


0j 9j 





2 i 














In Appendix A, the theorem is stated in a more general 
fashion along with the proof of it. In other words, allow¬ 
ing the distribution of phase cuts to deviate from the uni¬ 
form distribution results in non-classicality tests, which 
are equally efficient for the squeezed and unsqueezed sin¬ 
gle photon states. The introduction of the squeezing thus 
does not degrade the test and the results of Fig. 2 are 
also applicable to the squeezed single photon state. This 
is in contrast with the reconstruction method based on 
standard inverse Radon transformation, where squeezing 
would introduce additional ambiguities and/or a reduc¬ 
tion in the statistical significance of the results as dis¬ 
cussed below. 


III. INVERSE RADON TRANSFORM 
A. Introduction 


w 


R 



R W 


?\ </ Rt 

rIrip 


FIG. 5 (colour online). Action of the operators R and R'. 
The inversion cannot be achieved directly through the appli¬ 
cation of these two operators, i.e., there is a final step missing 
represented by the question mark. Since the latter is hard to 
achieve, the idea of the filtered back-projection is to rely on 
R 1 for making statements about W(x,p). 


Above, we have presented a simple and general non- 
classicality test that detects negativity of a phase space 
distribution W (x, p) resulting in the breakdown of any 
classical description of the system. This objective could, 
however, also be achieved through standard state recon¬ 
struction methods based on the inverse Radon transform. 
As already stated in the previous section, we cannot di¬ 
rectly measure the distribution W ( x , p) with arbitrary 
precision. What we can measure is the probability distri¬ 
bution of the quadrature operators Qg , which is given by 
the marginal distribution known as the Radon transform, 


RW(6»,s) 


W(x,p) S(s — x cos 6 — psind)dxdp. 


(25) 

Formally, RW(0,s) is an integral transform obtained 
from integration along the rays perpendicular to the cut 
characterized by 9 at distance s from the origin. Fig. 1(b) 
shows graphically how this transform works. The inverse 
process is referred to as the inverse Radon transforma¬ 
tion, which allows to extract the properties of the un¬ 
derlying distribution W(x,p ) from RW(0, s). We will 
now consider the use of the inverse Radon transform to 
show non-classicality. Furthermore, we will compare it 
to the simple test discussed above in order to determine 
advantages and disadvantages of the two methods for an 
experimental demonstration of non-classicality. 

The inversion of the Radon transform is not a simple 
task, but because of its importance in image reconstruc¬ 
tion, there is a wide range of techniques to accomplish 
it. For our purpose, the reconstruction of the complete 
distribution W(x,p) is not necessary, since we are only 
interested in the most efficient way to establish negativ¬ 
ity of W(x,p). This negativity of W(x,p) only occurs in 
limited regions of the phase plane and we thus do not 
need the full information. This motivates resorting to 
the so-called filtered back-projection theorem [20, 21], 
which provides a technique to evaluate the convolution 


of W(x,p) with another function on the phase space. 
Here, we will only present the main ideas of the theorem 
and refer to Ref. [20] for details. The Radon transform 
identifies a linear operator R on the functional space of 
well-behaved functions defined on the phase space, to the 
space of functions defined on the so-called unit cylinder, 
spanned by (0, s). For instance, W(x,p) belongs to the 
first functional space while RVF(0, s) belongs to the lat¬ 
ter. Both spaces can be equipped with suitable inner 
products and accordingly, if the Radon transform opera¬ 
tor is defined on functions like W(x,p), the definition of 
the adjoint operator Rf follows. It can be proven [20, 21] 
that the action of Rf on functions defined on the unit 
cylinder is given by 

, 1 f + % 

Wu>(x,p) = — / uj( 6,xcosO +psm6)d6. (26) 

nJ -f 

Despite this formula makes the adjoint Radon transfor¬ 
mation easy to compute, it does not solve the problem 
of determining the distribution W(x,p) because the ad¬ 
joint Rf is not the inverse of R (see Fig. 5). Instead of 
inversion, we can, however, use the adjoint to determine 
properties of W(x,p) from the filtered back-projection 
formula, which states that 



x, po — p)W(x,p)dxdp = 


RW (0, xq cos 9—po sin 0—s) u(9, s) dsdd , 

(27) 


provided that the filter, F(x,p ), and its kernel, ui(9,s), 
are related by the adjoint Radon operator, that is 

R f w (x,p) = F(x,p). (28) 

The advantage of the filtered back-projection formula 
in Eq. (27) in view of applications is straightforward. By 
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integrating the right hand side, where RIF (0, s) are the 
measured distributions and u>(9, s ) is a known function, 
one attains the convolution of the unknown W{x,p) and 
the known filter F(x,p). For state reconstruction, the 
idea is to avoid a direct cumbersome inversion of R by 
using Eq. (27) to get a coarse graining through an ad hoc 
designed filter. It is, however, important to realize that 
the filtered back-projection theorem shifts the problem of 
the inversion of the operator R to the problem of finding 
the kernel satisfying Eq. (28) for a given filter F(x,p). 
This problem, also called Cho’s problem in the literature, 
is highly non-trivial. For instance, a solution of Cho’s 
problem for F(x,p) = J^( x,p ) with & given in Eq. (3), 
is not known in general. 

For our purpose, the filter F(x,p ), will play the role of 
a test function and Eq. (27) gives us a recipe to compute 
its mean value, possibly translated in the plane, with re¬ 
spect to the unknown distribution W(x,p). Hence, sim¬ 
ilar to the simple test above, the non-existence of a pos¬ 
itive probability distribution W{x,p) can be proven by 
the right hand side yielding a negative value for a non¬ 
negative filter F(x,p). The difficulty of solving Cho’s 
problem, however, limits the choice of test functions that 
we can consider. 


B. Application 


In the previous non-classicality test, we looked for a 
violation of the condition in Eq. (3) for the squeezed and 
unsqueezed single photon state. Following a similar path 
of reasoning, we argue that a state is non-classical if it 
is possible to get a negative outcome by performing the 
convolution of W with a non-negative filter function F 
(see Eq. (27)) through the back-projection formula with 
a kernel oj as the solution to Cho’s problem. In the fol¬ 
lowing, we will refer to the filter F as the test function. 
A choice of filter is the normalized characteristic func¬ 
tion of a disc, for which Eq. (28) has been solved in the 
literature [22], 


Fa(x,p) 



(. x 2 +p 2 ) < a 2 
(, x 2 +p 2 ) > a 2 


(29) 


This filter is indeed non-negative and has a compact sup¬ 
port so that it can only have a negative expectation value 
if W(x,p) takes on negative values. This filter is not 
usually employed in tomography because it is not robust 
against noise. Other filters employed in reconstructions 
are more robust against noise, but they are not non¬ 
negative. The convolution with a negative filter would 
introduce an ambiguity in the interpretation of a possible 
negative outcome, since one would have to determine if it 
came from the distribution W or the filter F. We, there¬ 
fore, choose to consider a completely positive filter of the 
type in Eq. (29) despite its inherent instability. Ideally, it 
would be better to employ a more smooth non-negative 
filter, but to our knowledge no simple solutions to Cho’s 


problem has been determined for such filters. Using a 
circular disc is, however, only a good test function for 
the unsqueezed single photon state, which is also rota¬ 
tional invariant. The geometry of the squeezed state (see 
Fig. 4) suggests that it is better to use a stretched filter 
F a ,i(x,p) = F a (lx,p/l) for this state, provided that we 
can determine the corresponding kernel. Such a stretched 
filter will have a better overlap with the negative areas 
in the Wigner function of the squeezed state and we thus 
expect it to produce statistically more significant results. 
Note that for the sake of generality, we here allow for a 
stretching parameter l. which might be different from the 
squeezing parameter A. 

The solution to Cho’s problem for F a (x,p) given by 
Eq. (29) can be found in Ref. [22] and reads, 

«,„(»)= / -j (30) 

\/ 1— ( a / s ) 2 J ' ' 

In Appendix B, we derive a method to determine the 
kernel for any linear transformation of the filter. In par¬ 
ticular for the stretched filter F at i(x,p) we find the cor¬ 
responding kernel 

s) = v?&0) UJa (u&Fj) ’ ^ 

where we have defined the function 

u(9 , /) = C ° S ^ \/l + F tan 2 9. (32) 

In order to use the back-projection formula, we need 
to determine the integral on the rhs of Eq. (27). Since 
we know that the Wigner functions in consideration are 
negative around the origin, we set (xq,pq) = (0,0) and 
put l = A. The motivation for the latter is that if the filter 
covers the largest possible area in phase space, the error 
associated to the filtered back-projection is reduced, at 
least for rather small areas, since we are coarse graining 
over a larger area of the phase space. On the contrary, 
a filter limited to a small region around the minimum 
of W(x,p) will give a more negative result and thus a 
larger violation of classicality. The optimal filter is thus 
obtained as a trade off, which suggests a filter matching 
the symmetry of W(x,p). Consequently, we put l = A. 

We define ( F a> \) as being the lhs of Eq. (27), be., 

(Fa, A> = - / +2 / R W x (e,s)u a> x(e,3)d3d9. (33) 

nJ-f Jr 

Similar to Eq. (8), this equation expresses the mean 
value of the positive function F over the phase space 
in terms of measurable quantities since, for fixed 9 , the 
function RWa ( 9 , s) is just the measured distribution of 
the quadrature Qg. Nonetheless, the calculation of the 
mean value of the test function in Eq. (33) involves an 
integration over both the values of the quadrature and 
the phases. For the elementary test, only a finite set 
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of at least 21V + 1 phase cuts was required and a nega¬ 
tivity was already seen for IV = 4 in the single photon 
state. Here, on the other hand, evaluation of the filtered 
back-project ion over a certain finite number of angles will 
necessarily involve some ambiguities since it is a priori 
not clear how many angles are required to get a certain 
precision. Since we do not want to assume a priori in¬ 
formation about the phase space distribution, a proper 
analysis ought to include a study of the largest possible 
error introduced by restricting the number of angles. 

As opposed to estimating the error introduced by the 
finite number of cuts, we will in the following section 
focus on an alternative strategy. We will refer to this 
as a Monte Carlo numerical quadrature, where we also 
consider the measurement angle as a random variable. 
In this way, the number of cuts is equal to the number 
of measurements and we thereby gradually remove the 
ambiguities from the finite number of cuts as the number 
of measurements increases. We find that this strategy is 
conceptually simpler to understand, but an alternative 
and more natural strategy of using a fixed number of 
phase space cuts might be desirable experimentally. In 
the main text, we will only discuss this latter strategy 
briefly and leave most of the discussion to Appendix B 
where we show that it gives similar results although at 
the cost of introducing unpleasant ambiguities. 

For an experimental test, one is concerned with the sta¬ 
tistical error in the evaluation of Eq. (33). This depends 
on the specific strategy adopted but there is a general 
issue in all the test we consider, which is related to our 
choice of kernel to. The discontinuity of the kernel (see 
Eq. (31)) causes the variance of (F a , x ) to be infinite be¬ 
cause u) a '\(Q,s) is not square integrable for any A. The 
singularity is approached as e -1 for e = |s — «o I ~> 0, 
where so is the point where the singularity occurs. This 
infinity is simply a mathematical ambiguity and has no 
physical meaning. It can be avoided using a subsidiary 
kernel 


^a,X,e(9 > ^) 


0 

ui a ,x{0,s) 


a < |s|/u(A, 9) < a + e 
elsewhere, 


(34) 


which corresponds to removing a minimal interval in the 
values of the quadrature close to the singularity. This 
results in a small systematic error in the estimation of 
(■ F a ,x ), equal to 


A a (e) = 


/*+■?£■ p (a+e)it(A,0) 

' — Ja u(X,6) 

pa-\-e 


RWa ( o , s) 


|w aiA (6 , ,s)| d9ds 


pa- 

/ RRA=i (t) |w Q , A =i(t)| dt. 
J a 


(35) 


Note that this quantity cannot be directly computed from 
experimental data. Nonetheless, an upper limit on this 
error can be obtained empirically, by replacing Rif with 
its maximum measured value. This leads to an error of 
the order O ^ . The variance for the estimate of (F a: \) 
is thus finite, but proportional to | log e|. In the following, 


the use of the the subsidiary kernel w aj A,e(0, s) will always 
be intended even though we will not explicitely write the 
subscript e. 


C. Monte Carlo numerical quadrature 


We now introduce the Monte Carlo numerical quadra¬ 
ture for the double integral appearing in the back- 
projection formula. The goal is to divide the 
integrand into two parts, RW x (9,s)uJa, X (9,s)/n = 
p x (9, s) I a , X (9, s), where p x (9,s) is considered to be a 
probability distribution for two random variables 9 and 
s, and I a ^\(9,s) is a suitable integrand reminder. Note 
that in the experiment, only the variable s, which charac¬ 
terises the quadrature value, is a truly random variable 
resulting from the probabilistic nature of the measure¬ 
ment, whereas the phase angle can be controlled exper¬ 
imentally. Nonetheless, we assume that the phase angle 
is also varied randomly so that, for each shot of the ex¬ 
periment, the experimentalist picks the phase angle ac¬ 
cording some probability distribution C(9). According to 
Eq. (33), we can then express the filtered back projection 
as 

{F a ,x) = [ [ I a ,x{0, s )p\{9,s)d s d0. (36) 

J - f Jr 

with 

f Px (9 lS ) = RW x (9,s)C(9) 

) I ,(Q s )= V 67 ) 

where C ( 9 ) is any distribution of cuts normalized to unity, 
i.e., f C(9)d9 = 1. At this point, all distributions of cuts 
are logically equivalent and will give the same result for 
{F a , X }, but as we discuss below the choice of distribution 
will make a difference for the efficiency of the numerical 
quadrature. 

In an experiment, we assume to obtain a set of M 
outcomes of the form {(0j,Sj)}, distributed according to 
p x (9,s), so that the average of the test function can be 
estimated as 

1 M 

(Fa, a) — y s i)- (38) 

i— 1 


The estimate can then be used to check the classicality 
condition in Eq. (3) . For the statistical significance of the 
test we, however, need to include the error in the Monte 
Carlo quadrature. The variance of the estimate (F a , x ) 
can be expressed as 

, m 2 

4„ lX = JmZTY) E - <^A>) (39) 

Accordingly, the relevant quantity to consider for the vi¬ 
olation is the ratio 

(Fg,x) 

&F a , x + A a(e)’ 


( 40 ) 






11 


which accounts both for the statistical uncertainty er p a x 
and the systematic error A a (e). Since we now consider 
not only the mean but also the statistical uncertainty, 
such a ratio depends on the distribution of cuts C(9), 
which plays an essential role in the minimisation of <Jp a x - 
This is discussed below and in Appendix B, where we 
consider both the squeezed and unsqueezed single photon 
state. Now there is no longer a sensible way to remove 
the dependence on the number of measurements M from 
this ratio. In fact, the asymptotic behaviour for large M 
is no longer ~ \J~M as for the simple test because there 
is a the trade off between <jp a A , diverging as | log e|, and 
A a (e) decreasing as e 2 . Therefore it is not convenient to 
rescale it by yM and take the limit of M —> oo as we 
did for the simple test (see Eq. (14)). In order to achieve 
a comparison with the simple test discussed above we 
consider the ratio 


K(M) = 


{&) + A a (e) 

& & (F a ,x) 


(41) 


which describes the ratio between the statistical certainty 
of the tests, i.e., the ratio of the number of standard 
deviations with which the tests violate the classicality 
criterion. A value 7 Z(M) > 1 (7 Z(M) < 1) thus signifies 
that the simple test (inverse Radon method) has a larger 
statistical certainty. 



FIG. 6 (colour online). Plot of the ratio 1Z(M) of the statisti¬ 
cal significance of the elementary test and the inverse Radon 
method for the single photon state. The comparison is with 
respect to the best elementary test for TV = 16, m = 17. The 
value 1Z(M) > 1 indicates that the elementary test is statisti¬ 
cally more reliable. The almost linear part of the graph shows 
a nearly logarithmic growth for high M; the low M region de¬ 
viates from this behaviour due to the y/M — m dependence 
in <tjt. The dashed vertical line indicate the threshold for the 
inverse Radon method to violate the non-classicality by two 
standard deviations. 


D. Single photon state. 

First, we consider the back-projection formula for the 
unsqueezed single photon state. The rotational invari¬ 
ance of the Wigner function W(x,p) motivates a phase 
cut distribution of the form 

m = (42) 

7r 

which corresponds to a uniform, flat distribution of the 
phase cuts. As shown in Appendix B, this is the optimal 
distribution of cuts for minimizing ap a A . Consequently, 
we have 1(9, s) = Tru! a ,\(0,s). 

For given M, we have minimized the value of 
<j f ^+a (e) w ith respect to the remaining free param¬ 
eters a and e, keeping A = 1 because of the rotational 
symmetry. In the minimisation, we numerically com¬ 
pute the double integral that estimate the value of the 
mean (F at \) and variance crF aX for given a and e and 
analytically compute the upper bound of the systematic 
error A a (e). Consequently, we can evaluate the ratio in 
Eq. (40) and run a numerical optimisation to minimize 
it. We compare the result of the minimisation to the 
best elementary test we could achieve with TV = 16 (TV 
is the polynomial degree) by calculating the ratio 1Z(M). 
The results of this is shown in Fig. 6 as a function of M. 
From the optimisation we find that the systematic error 
decreases as a power law and this causes the variance 
of the kernel to grow logarithmically. Consequently, the 


quantity TZ(M) does not settle to an asymptotic value, 
but keeps increasing as a function of M because of the 
trade off between the systematic error and variance of the 
kernel. For all values of M we find R > 1 which means 
that the elementary test is statistically more significant. 
The inverse Radon is closer to the elementary test for a 
small number of measurements. Note that the very small 
M region cannot be investigated since crjr is defined only 
for M > m, the least number of cuts being to = 17 for 
TV = 16. This region is, however, not very interesting 
since it lies in the region where the inverse Radon is not 
capable of producing a violation by two standard devia¬ 
tions (dashed vertical line in Fig. 6). 

E. Squeezed state. 

We now consider the squeezed single photon state, 
where an explicit optimization in the choice of C(9) is 
crucial for obtaining an effective test. As shown in Ap¬ 
pendix B, the naive approach of choosing C(9) as given 
by Eq. (42) leads to a poor scaling of the variance as ~ A 2 
(~ A -2 ) for A 3> 1 (A -C 1). This result shows that for 
the squeezed state a uniform distribution of phase cuts is 
not desirable since the squeezed state has a highly non- 
uniform angular distribution. It is thus better to have 
more measurements in regions of phase space with a rapid 
variation. To find a more optimal distribution C(9), we 
use the following result from the theory of Radon trans- 
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Hence the same conclusions drawn about the statistical 
significance of the elementary and inverse Radon tests for 
the single photon state also applies to the squeezed single 
photon state. 


F. Finite Number of Cuts 

In experiments, a strategy more widely used than the 
Monte Carlo quadrature is to have a small set of cuts 
but many measurements for each cut. In order to study 
this situation, it is more convenient to work out the right 
hand side of Eq. (27) in a different way than for the Monte 
Carlo quadrature. This is done in two steps. First, we 
introduce the function 


FIG. 7 (colour online). Logarithmic plot of x) as a func¬ 
tion of 9 for some values of the squeezing parameter A. The 
region that contributes the most to the integral in Eq. (44) 
gets narrower as lambda increases. 

form [20] 


g a ,x{0)= f RW\ (0,s)u a ,\(0,s)ds, (46) 

Jr 

which is assumed to be accurately estimated by having 
many measurements per cut. The integral over the phase 
cuts is then written as a summation over a convenient 
choice of cuts {9j} with weights Wj 


RVCa (9, s ) 


u(A, 9) 


HW 


it(A, 9) 


(43) 


together with the transformed kernel for the stretched 
filter given in Eq. (31). Using these two results, we find 
that 


(Fa,\}= f +2 I RW A (M h,a ' A(g,s) dsde 


1 


! (M) 


RWa=i (t) dtd9. 


(44) 


Eq. (44) shows that not all angles contribute equally to 
the final result. Because of the use of a stretched filter 
most of the weight is focused on a small subset of angles. 
This is illustrated in Fig. 7 where we plot u2 ^ e as a 
function of the phase 9 for various values of A > 1. As 
shown in the figure, u2 ^ e ^ gives more weight to phases 
close to 9 = 0 {9 = tt) for A > 1 (A < 1). Since it is im¬ 
portant to get an accurate sampling of the most relevant 
part of the integral, it will be advantageous to choose a 
probability distribution for the phase cuts C(9) which is 
higher around this region. As shown in Appendix B, the 
optimal distribution of cuts is indeed 


m = 


i 

7tu 2 (A, 9) ’ 


(45) 


which coincides with Eq. (42) for A = 1. Furthermore, we 
show in the appendix that the relevant ratio (e) 

does not depend on A for this choice of phase cut dis¬ 
tribution, i.e., it is thus equal to the one found for the 
rotationally symmetric state. This is completely equiv¬ 
alent to the analogous result found for the simple test. 


9a,\{9)d9 ->• 5>, flo , A (A,). (47) 

{<hl 

The discretisation of the integral above, formally called 
numerical quadrature, is the crucial step in a strategy 
with a finite number of cuts. This choice, however, gives 
rise to ambiguities because we cannot completely know 
the function g a \ ( 9 ) since it is sampled with only a few 
phase arguments {9j}. The first ambiguity is in the 
choice of the number of phase cuts. We cannot a priori 
determine a minimum number of phase cuts required as 
in the elementary test. Regardless of the number of cuts, 
the full expression on the left hand side of Eq. (47) in¬ 
volves an integral. Hence one can always wonder if adding 
more phase cuts could reveal narrow features missed by 
the finite number of cuts. It is therefore in principle 
necessary to constrain the influence of such narrow fea¬ 
tures on the final result to convincingly demonstrate non- 
classicality. Another ambiguity is in the estimation of the 
error due to the numerical quadrature, which cannot be 
easily determined unless some other information about 
9a,\ [9) is known. 

Despite of these ambiguities, it is still reasonable to ask 
whether we can determine an efficient numerical quadra¬ 
ture that uses a low number of phase cuts compared to 
the total number of measurements. In other words, we 
want to understand whether the inverse Radon method 
can provide a way to choose a suitable small set of phase 
cuts if some further sensible assumptions are made. To 
give an answer to this question, in the analysis carried 
out in Appendix B, we assume ab initio that the states 
considered are the single photon state, described by a 
symmetric W\[x,p ), as well as the squeezed single pho¬ 
ton state. For our theoretical analysis, such an hypothesis 
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relies ultimately on the Wigner function to be the cor¬ 
rect description of the state, which is somewhat logically 
problematic since this information is exploited to carry 
out the method, not just to optimise it. 

We show in Appendix B that this strategy can also be 
carried out in such a way that the statistical strength of 
the test is the same for both rotationally symmetric and 
squeezed states. As the squeezing increases, however, 
more cuts are needed in order keep the uncertainty of 
the test below a certain threshold and these cuts should 
be highly asymmetrically distributed. Choosing such an 
asymmetric distribution can be justified by the asymme¬ 
try of the filter function, which is highly peaked around 
a narrow distribution of angles. On the other hand, it is 
hard to quantitatively account for the absence of further 
errors due to the unknown angular dependence of the un¬ 
derlying physical state, and this introduces an ambiguity 
in the test. 

In essence, closing such “loop-holes” in the argument is 
what is achieved with the elementary test, which is ex¬ 
plicitly designed to deal with a finite number of cuts. Due 
to the ambiguities discussed here, making a more direct 
comparison between the elementary test and a strategy 
with a finite number of cuts is highly non-trivial and is 
thus beyond the scope of this article. 

IV. CONCLUSION AND DISCUSSION 

In conclusion, we have generalised the non-classicality 
test presented in Ref. [6] and developed a framework for 
finding the optimal test function for a given quantum 
state. In the generalisation of the test, we consider test 
functions which are arbitrary polynomials in the phase 
space. Thereby, we allow test functions which do not 
have rotational symmetry and are therefore more suited 
for testing the non-classicality of quantum states which 
do not possess this symmetry. As a specific example, we 
have found the optimal test function for a squeezed single 
photon state and shown that it results from a simple 
transformation of the rotationally invariant test function 
considered in Ref. [6]. 

Furthermore, we have performed a comparison be¬ 
tween the simple non-classicality test and more conven¬ 
tional tests based on the inverse Radon transform. In 
order for the comparison to be fair, we have extended 
the strategy of inverse Radon transformation based on 
the filtered back-projection formula so that it is more 
optimised for squeezed states. The constraints of the fil¬ 
tered back-projection formula and the limited number of 
known completely positive filter functions meant, how¬ 
ever, that only a rather restricted class of test functions 
could be investigated (as compared with the elementary 
test, where there is a lot of freedom). In addition, the 
filtered back-projection formula needs to be handled with 
care because it introduces problems with divergencies 
and ambiguities from introducing a finite number of cuts. 
With a Monte Carlo quadrature method the ambiguities 


from the finite number of cuts can be avoided, but this 
requires the use of as many quadrature cuts as the num¬ 
ber of measurements. Reducing the number of cuts in¬ 
troduces undesired ambiguities, which may question the 
whole outcome of the test. In particular, these ambigui¬ 
ties are severe for strongly squeezed states where it is nec¬ 
essary to have a large number of cuts in a small angular 
interval where the integrand varies rapidly. Without any 
prior information about the state, such a choice of dis¬ 
tribution introduces ambiguities, which make it hard to 
give a convincing proof of the non-classicality of a state. 
All of these problems are avoided with the elementary 
test. In addition to this, we find that the elementary 
test is statistically stronger than the methods based on 
the inverse Radon transform (at least compared to the 
version based on the Monte Carlo quadrature method). 

We thus find that the simple test is more suited and 
straightforward than the inverse Radon transform for ex¬ 
perimental demonstration of non-classicality. The ad¬ 
vantage of the simple test originates from the fact that 
it is naturally designed as a non-classicality test under 
the constraints imposed in a typical experiment. In par¬ 
ticular, the elementary test is specifically designed to 
work with a finite number of quadratures being mea¬ 
sured. On the contrary, the inverse Radon transform is a 
more general method, which allows to gain more informa¬ 
tion about the physical state but it is less suited for this 
specific situation. The statistical weakness of the latter 
could, however, be due to the fact that the kernels con¬ 
sidered here are not ideal due to their discontinuity. On 
the other hand, the smooth filters considered in the liter¬ 
ature that grant a better convergence are not necessarily 
completely positive. A negativity observed with such a 
filter would therefore inevitably introduce the question 
of whether the result was caused by the negativity of the 
filter and this could compromise the test. 

We have considered the simple test for a quadrature 
squeezed single photon state but the general framework 
developed here can be applied to any quantum state or 
any physical system. In particular, it could be inter¬ 
sting to use the test to demonstrate non-classicality of 
opto-mechanical systems [10, 1 ]. Currently much ef¬ 
fort is devoted to reaching the quantum regime for these 
systems. To convincingly show that such systems are 
outside the realm of classical physics, it is essential to 
demonstrate a contradiction with classical physics along 
the lines discussed here. Also it could be interesting 
to apply the method to investigate, the negativity of 
macroscopic quantum states such as large Schrodinger 
cat states and potentially use this as a measure of non- 
classicality in such systems. Finally, one could imagine 
using the elementary test method as a general tool for 
tomography without resorting to the heavy machinery of 
inverse Radon transformations. Since we find the method 
to be more reliable and less ambiguous, one could imagine 
other situations where this method could be a desirable 
strategy for tomography. 
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Appendix A: Test function method 

In this appendix, we present the details of the results 
on the simple test method discussed in the main text. 
First, we prove, in three intermediate steps, a slightly 
more general version of Theorem 1. The first step is to 
establish that for a rotationally symmetric distribution 
the average value of any test function can be expressed 
as the average value of a radial test function. The second 
step is to identify the tests that minimise the variance 
for a given value of the mean. Finally, the explicit con¬ 
struction of a radial test which minimises the variance is 
provided. The proof of the theorem is then obtained by 
combining these three partial results. Second, we prove 
Theorem 2 concerning squeezed states. 

1. General formulation. 

We will take into account a slightly more general form 
of the test function than in the main text in order to 
derive the proofs of the theorems. To begin, we consider 
the test function to have the form 

2N 

^(x,p) = l+^2Pi(x,p) (Al) 

i=1 

where Pi(x,p ) is a homogeneous polynomial of degree i 
defined on the phase space of x and p. Moreover, all the 
test functions of the form in Eq. (Al) are implicitly taken 
to be non-negative, thus we are more general than in the 
main text where we only considered square polynomials. 
Note that the polynomials Pi can be expressed in terms of 
the functions x l ~ k p k to recover the notation employed in 
the main text of the paper. By assumption, we expect & 
to be in the algebra of observables, formally defined as the 
commutative algebra generated by all the quadratures. 
In particular, for fixed i, any Pi(x,p) belongs to a linear 
subspace of this algebra, together with Q l s for any 9. Any 
set of at least (z + 1) functions of the form Q l g for different 
9 is a (possibly overcomplete) set in this linear subspace. 
Accordingly, given a set of j = 1,.., m different cuts {9 : j}. 
with m > 2N + 1, all the polynomials Pi(x,p) involved 
in Eq. (Al) can be expressed as 

2AT+1 

Pi(x,p) = ^ Ti.jQi., (A2) 

j=i 


with Ti-j being real-valued but not uniquely defined. We 
will express a test, T as 

IT = (N, Pi, 9j, Ti-j) , (A3) 

where the non-negativity of the test function is intended, 
but not explicitly stated. This is the most general class 
of tests we consider and we refer to this class as . A 
test of the form in Eq. (A3) belongs to this class, Te?, 
if & as defined in Eq. (Al) is non-negative and Eq. (A2) 
holds. 

The rotated coordinates, Q$ and P^ are defined as 



and correspond to the quadrature Q^ and the rotated 
canonical momentum P^. Accordingly, we define the ro¬ 
tation of polynomials by 

Pi[4>\{x,p) = Pi{Q<j,,P<j>)- (A5) 

We refer to a test as a radial test and denote it by T ra d, 
if Pi[cj>\(x,p) = Pi{x,p) \/cj) £ R, for all i = 1, ..,2 N + 1. 
Consequently, in this case, the test function, ^, is only 
a function of the radius. The class of such radial tests 
is denoted Ifrad and it is obviously a proper subclass of 
^. To complete the set of notations, it is sensible to 
explicitly express the dependence of (&) / o& in terms of 
the free parameters of the test. Therefore we shall write 

( ^) 

^-=g[T,M], (A6) 

<7J? 

so that the actual figure of merit Q is simply recovered 
as e[T] = liniM^oc 

2. Rotationally invariant states. 

As hinted above, even though our objective is the opti¬ 
mal <?[T, M] for any M, we do not yet consider the direct 
optimisation with respect to the free parameters of the 
test. Instead, we first want to understand the structure 
of the class ^ and its relation to ^ ra d, for a rotationally 
symmetric W(x,p). The first result concerns the mean 
value (&) and, in particular how the class ‘rf and ffrad 
are related under the condition that all (P t ) are given. 

Proposition A.l. Assume that W[x,p) is rotationally 
symmetric and consider a generic test T £ to such that 
AP = 1 + Y^i=i Pi( x ->P) the test function. There always 
exists a radial test T ra d £ c 6.’ r ad with test function of the 
form &rad = 1 + YH=\Pi( x iP) L where Pi(x,p) is rota¬ 
tionally symmetric and fulfils (Pi) = (Pf) for any i, so 
that (TP) = (^ r ad) ■ Furthermore, all the radial tests that 
satisfy this condition differ from each other only in the 
parameters {9j} , {Ti-j}. 
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Proof. The proof is by construction. The given test is 
characterized by T = (N, Pi, 9j, T^j), and we define 

1 f 2n 

Pi = — / Pi [0] dd Mi = 1,2 N. (A7) 

27r J o 

Note that Pi is rotationally invariant since 


It can then be seen that a 2 = (j#’ 2 ) — (MMj) 2 = 

\ Y^i'~\ so that a final expression for a 2 ^ is 

achieved 



m 2N 


EE 

1=1 7. ?/ = 


9 


a 1 Tj-jTj'-j 
Mj - 1 ’ 


(A12) 


Pi [0 = ^ Pi [9 + <t>\ de = P % . (A8) 

Moreover, because of the rotational symmetry of W(x, p ), 
a simple change of variables to the rotated coordinates 
shows that (Pi [0]) = (Pi) for any 9. Thus, 

1 r 27r 

M‘) = 2 iJ 0 

= <a>-f r M={P .y (A9) 

Considering A? ra d = 1 +^i=i Pi( x ,p), it can be seen that 
the rotationally invariant test function is non-negative 
since 

2iV ^ «27r 

&rad(x,p) = 1 + ^ — / Pi [9] (x,p)d0 

i—1 ^ •'° 

= 2tt/ &[9\(x,p)d9, (A10) 

i.e., for fixed ( x,p ), it is the integral of a non-negative 
function AP ra d [$] ( x,p ) wrt. the variable 9. The statement 
of Proposition A.l follows immediately. Note that no 
consideration of 9j or T i;j was necessary because they 
are irrelevant for the mean value of &. 

□ 

Above, we have proven that for a rotationally sym¬ 
metric distribution, W(x,p), and a non-rotational test, 
we can always find a rotationally symmetric version of 
the test, which gives the same average. We will now 
show that the rotationally invariant test will also have a 
smaller uncertainty in the mean value. To do so, we ini¬ 
tially study a more general problem, neglecting for now 
the non-negative property of the test functions. Since we 
consider a more general class of functions, we will at first 
(see Proposition A. 2) achieve results without considering 
the requirements for proper test functions, and thereby 
of no direct physical significance. Later we will, however, 
show (see Proposition A. 3) that it is indeed possible to 
find a a physically meaningful test fulfilling the condi¬ 
tions of a optimum in the broader class, and thus also an 
optimum in the restricted class. 

First, we give a suitable expression for cr^ that ac¬ 
counts for the rotational symmetry of W(x,p). Since 
(' Q l g ) = ( Q l ) does not depend on 9 for any i, we shall 
introduce the covariance matrix of the quadratures, 

Er' = (q^'} - (o i ) (</) (ah) 


where {Mj} is the number of measurements for the j’th 
cut and ]Tb Mj = M. 

We now perform an optimisation in the class of func¬ 
tions of the form & = 1 + J2i jTi;jQg., with the only 
constraint ■ Tj;,- = Ki. This constraint on Tj-j origi¬ 
nates from the identity in Eq. (A2). In the case of rota¬ 
tional symmetry, where ^ Q l g ^ Q ^ does not depend 
on 9j, this implies that 

m /p\ 

= W> = Kh (A13) 

where K, are constants that depend only on the values 
(Pi) and no other free parameters of the test. By opti¬ 
mising the variance subject to this constraint, we thus 
optimise over the class of tests giving a fixed value of the 
mean value. Since the true optimum must also be the 
optimum for the particular mean value obtained, proper¬ 
ties proved under this constraint also applies to the true 
optimum. The optimisation is done keeping N, M and m 
fixed, since these parameters do not depend on whether 
the test is radial or general. We can thus state and prove 
the following result. 

Proposition A.2. Assume AP to be any generic poly¬ 
nomial function of order 2 N over an arbitrary set of m 
quadratures, that is 


m 2N 


^ = i+EE^A 


(A14) 


i=1 i=i 


We assume that all the quadratures Qg j are independent, 
identically distributed random variables and that T^j are 
only constrained by the relation Y^jLi 7-,j = dt fol¬ 
lows that the minimum value of Oj?, after optimising Ti-j 
and keeping fixed N, M and m, is 


2 1 Ei iLi d^Ki-Ki’ 

2 


M- 


(A15) 


and is thus independent of the particular choice of {Mj} 
satisfying y ■ Mj = M . Furthermore, if {Mj} are all the 
same, the bestTi-j is independent of j, i.e., T^j — Ki/m, 
m being the number of cuts. 

Proof. Consider a distribution of measurements {Mj}. 
For simplicity, we define Mj = Mj — 1 so that JT Mj = 
M — m. The constrained optimisation problem is then 
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studied using the Lagrangian function 

2N m 

i=1 J=1 



which only depends parametrically on {Mj}. The opti¬ 
mal solution is determined by 


the optimum in Proposition A.2. To this end, it is sen¬ 
sible to fix Mj to be all the same, so that the statistics 
per cut is unbiased. Then we will show that we can con¬ 
struct radial polynomials according to Eq. (A2) in such 
a way that Tjj is independent of j. Note that the only 
radial polynomial we can construct must be functions of 
(. x 2 + p 2 ) to retain rotational symmetry. 

Proposition A.3. Assume T^j to be independent of j. 
Then the polynomials Pi given by the identity in Eq. (A2) 
are radial polynomials provided that the phase cuts are 
uniformly distributed. In particular, the identity 


{ TEAL — \^ 2N 3"'T'; j \ _ n w- ■ 

dTi-j ~ E»'=l Mj Al ~ U V *’ 3 

ET=i T fJ-Ki=0 Vi ( Al7 ) 

EjLi Mj=M~ TO. 

It is seen from the first equation that the optimal solution 
is such that T i; j/Mj does not depend on j, i.e., it is 
independent of the choice of phase cuts. From this, we 
find that the solution fulfilling all three conditions is 


T ■ = 


Kj Mj 
M — TO 


(A18) 


(x 2 +p 2 ) n =T 2n Y / Q% (A20) 

3 =1 


holds, where T 2 n 
n + 1. 


2 n 
n 


-l 

— and 9i = — with to > 

m J m — 


Ti 


Note that we have dropped the j index in the coefficient 
= T 2 n and we have used the index n in place 




of i to avoid confusion in the subsequent proof of the 
statement, where the imaginary unit is used. 


Then, the constrained stationary value of the function in 
Eq. (A16) is 


Proof. Since Tj-j do not depend on j , we consider the 
following sum 



2N <7* 

i,i '=1 y 




M — TO 


(A19) 


It can be proven that it is a minimum by noting that the 
Hessian of Eq. (A15) is the tensor product of the covari¬ 
ance matrix and the diagonal matrix containing 1 /Mj on 
the diagonal. Since both of these matrices are positive- 
semidefinite by definition the found extremum is a mini¬ 
mum. Finally, it follows from Eq. (A18) that when {Mj} 

are all the same, the ratio is exactly 1 /to. 

□ 


The crucial part of the theorem is the independence 
of the covariance matrix in Eq. (All) on the choice of 
phase cuts, which is a direct consequence of the rota¬ 
tional symmetry of W(x,p). It follows that the choice of 
the cuts is irrelevant for Proposition A.2. Furthermore, 
the minimum of Cj? does not depend on the choice of 
{Mj}, provided that M is fixed. This is in contradiction 
with the physical intuition that it is not a suitable strat¬ 
egy to perforin most of the measurements on a single 
quadrature whereas only a few measurements are per¬ 
formed on the remaining quadratures. Note, however, 
that this result above is obtained in the broader class of 
functions only subject to the constraint T VJ = K,. 
Hence, for arbitrary choice of {Mj}, we might not be 
able to construct radial polynomials fulfilling the gen¬ 
eral requirement Eq. (A2) with the optimal given by 
Proposition A.2. As we will show below, it is, however, 
possible to find a proper test fulfilling both Eq. (A2) and 


£<$* (A21) 

3=1 

where for the sake of generality, the sum is extended over 
to cuts, with the suitable m to be determined. Using the 
exponential representation of the sine and cosine func¬ 
tions, we get 

Q% = ^ ( e ^ (* - ip) + e-i * J (* + *P)) 2n ( A22 ) 

1=0 ' ' 

so that 


3=1 


QS? = 



ip) 2n ~ l (x + ip) 1 . 

(A23) 


By choosing uniformly distributed cuts, we can turn the 
summation over j into a simple geometric series. Letting 
6j = the geometric series is 

m 

E (AM) 

3=1 
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and the exponential is different from unity if is not 
an integer. Accordingly, if to > n + 1, we have that 

m 

J2e i20 ^ n ~ l ) =m6i >n (A25) 

3=1 

Inserting this into Eq. (A23) gives 



and hence the statement of the proposition is proven. □ 

Note that only n +1 cuts are needed to express a poly¬ 
nomial of degree 2 n such as (a; 2 + p 2 ) . This is because 
odd powers of x and p are not involved. Therefore, for a 
radial test function of degree 2TV (corresponding to the 
square of a polynomium of degree TV), we only need a 
minimum number of TV + 1 cuts. 

We now have all the intermediate results necessary to 
prove the statement of Theorem A.l, which constitues 
the main result concerning the optimal Q[ T, M\. 

Theorem A.l. Assume that W(x,p) is rotationally 
symmetric and takes on negative values; then the opti¬ 
mal test, i.e., the test T;, est that minimises Q(T,M\ for 
any M belongs to the subclass ’ rad and the cuts are uni¬ 
formly distributed on the plane (see Eq. (A20)). 


Therefore, defining T rad = (TV, Pi, Oj, T v j) and using 
that = {^radj < 0, we arrive at 

Q[T rad ,M\<G[Y,M\. (A29) 

□ 

Thus, we have proven that whenever W(x,p) is rota¬ 
tionally symmetric, a rotationally invariant test is opti¬ 
mal. It should be noted, however, that the test described 
in the main text is not as general as the one considered 
in this appendix because the former is always a square 
polynomial to ensure non-negativity in a simple manner. 
Moreover, the map we have used 

-► Pi (A30) 

does not preserve the property of being a square polyno¬ 
mial for &, that is & rad = 1 + J2i Pi is n °t generally a 
square polynomial even when & = 1+JT Pi i s - Nonethe¬ 
less, the square root of AP rad can be efficiently approxi¬ 
mated by a high order polynomial in r 2 in the relevant 
region where it is peaked at the negativity of the Wigner 
function. This, together with the convergence shown in 
Fig. 2, motivates the choice of polynomials in r 2 . 

3. Squeezed states. 


Proof. Let us fix M and pick up a generic test T £ ^ 
characterized by (TV, P i; Oj, T^j), for which Q[T,M\ is 
negative. We will construct a radial test T rad £ ^’ ra d for 
which g[T rad ,M] < Q[T,M\. 

By Proposition A.l, we can find a set of radial poly¬ 
nomials Pi such that (Pi) = (Pi) and clearly ^ = 

^ rad ). By symmetry considerations, we know that 
only radial polynomials with even degree will contribute 
to the mean value. We now put Ify = T 2 n = ^Q 2 L ) ■ By 
Proposition A.3, we pick up a set of uniformly distributed 
cuts, Oj = with m being, for simplicity, the number 
of cuts in T, so that the identity in Eq. (A2) holds, i.e., 
P 2n = T 2n Q“§ n - Furthermore, it is seen that both T i;j - 

and Ti-j satisfy the same constraint 


We will now consider the behaviour of the optimal test 
for the more general squeezed state. 

Theorem A. 2. (Invariance under squeezing) Consider 
the best test To = (TV, Pi, Oj, Tj.j) £ ^ for the joint 
distribution W(x,p). Then the best test for the cor¬ 
responding squeezed state W\(x,p) = W(\x,p/\), X be¬ 
ing the squeezing parameter, is characterised by T\ = 
(tv, if>, ef\ 2g>) where 

Pr X \x,p) = P l (Xx,p/X), (A31) 


tan 0^ = (A32) 




3=1 



3 = 1 


{Pin) 

(Q 2n ) ‘ 


(A27) 


Consequently, both & = 1 + jPi jQg, and &rad = 
1 +fP, j Ti.jQ l g satisfy the hypothesis of Proposition A.2 
with the same Ify = It follows from the propo¬ 

sition that Ujir is the minimum attainable variance 
within that class or equivalently 

2 / 2 
°> ro d — 


rp(X) _ rp 

± i-,j — J -ro 


cos Oj 
cos Oj^ 


(A33) 


Proof. The proof is by construction of a bijective cor¬ 
respondence between tests for the unsqueezed case and 
tests for the squeezed one. First, we consider the repre¬ 
sentation of the new set of polynomials in terms of the 


(A28) 
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new set of cuts 

P t W (x,p) = P z (Xx,p/X) (A34) 

= ^ Ti-j ^cos 9j Xx + sin 9j ^0 


= E ^ ( cos6 '? A ) 1 +—^p 


In the domain (—7r,+7r), we can find new cuts 9^ ac¬ 
cording to Eq. (A32) and obtain 

P[ X) (x,p) = 'Y^,T i . j (a C0S ^ ) > ) (cos 0^ A) x + sin^ A) p) 


cos 9^ 


3 \ ~ 3 

'0\ 

3 


(A35) 


where is defined in Eq. (A33). 

We will now show that ^ = ^^(*0 and that 

g 2 ^ = g^ w where & is the test function for 

the unsqueezed (squeezed) state. We recall the result 
for the Radon transform in Eq. (43) stated below for 
convenience: 


RW a (9 , s ) = 


-R W 


(X,9) \u(X,9) 


(A36) 


where A is the squeezing parameter and it(A, 9) is defined 
as u(X,9) = ^^a/I + A 4 tan 2 9. It follows that 

(<70 = f + dss l KW x (9,s) 

' ' A J-CX) 

ds 


' — OO 
f» + 00 


u{X, 9) 

C+oo 


HIT 


v u(A,0) 

/ +oo 

dtf RW (0,t) 

-oo 


so that (j0 = («^ A 0 ■ 

We now wish to establish a similar result for the vari¬ 
ances, i.e., o'lr(A) = We have that 

1 2N 

= 2 E E 9xVj) ( A4 °) 

3 i,i '=1 

where the covariance matrix \g l x (%) = 

{J\Q l 8~j l ) — (Qh ) ) ^ as 4)6611 introduced. 

Eq. (A37) implies that g x ( 9j ) = u l+l (A ,9j)g vl so that 

1 2N 

= 2 E E = a%, (A41) 

j i,i' = 1 

because u l 0, dj A 0 = ^A cos Jj,) j . Note that the vari¬ 
ance for each cut a 2 is preserved. 

To conclude, we have explicitely constructed a corre- 
spondance between any test used in the case of no squeez¬ 
ing To and a test to use in the squeezed case T> such that 


g[T x ,M\=G[T,M\. (A42) 

This result establishes that the best test for the 
squeezed state is at least as good as the one for the 
unsqueezed state. The transformation applied, however, 
also works in the opposite direction, so that a more op¬ 
timal test for the squeezed state, could be unsqueezed 
and act as a more optimal test for the unsqueezed state. 
Hence there is a one to one correspondence between 
the optimal tests between the squeezed and unsqueezed 
states and the statement of Theorem A.2 follows. □ 


Appendix B: Inverse Radon method 


= u\X ,0)(q 1 }, (A37) 

where we have changed from the variable s to t = ^ . 

For the mean values, we find 



Since tan 00 = tA "/* , it is seen that u*^A, #j A 0 = 
(a and hence 

(^Ve t ^) £ (4 (A39) 

3 


In this appendix, we discuss the results of the filtered 
back-projection algorithm that have been quoted in the 
main text without proof. We generalise and prove the re¬ 
sult used for the linear transformed kernel (see Eq. (31)), 
and provide the proof that the distribution of cuts de¬ 
scribed for the Monte Carlo quadrature method is the 
optimal distribution. Finally, we discuss the alternative 
procedure with a finite number of cuts in detail. 

The result for the linear transformation of the kernel 
can be derived for the d-dimensional case. The planar 
case (d = 2) that has been presented in the main text will 
then be retrieved as a special case. The whole theory of 
the Radon transform can be formulated within where 
the Radon transform is achieved by integration over all 
hyperplanes in rather than straight lines in the plane. 
For a detailed but not too abstract discussion see for 
instance Ref. [20]. Here, we only focus on the adjoint 
Radon transformation. Let x denote a vector in and 
(• | ■) denote the standard scalar product. A non singular 
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linear transformation from R d into itself is denoted by L 
and we introduce the (d — l)-dimensional sphere, S d ~ x , 
as the set unit vectors in R d . The measure of the total 
set 5' i_1 is 


\S d ~ 1 \= [ (in. (Bl) 

J lies ' 1 - 1 

We now give the explicit functional form of the adjoint 
Radon operator RJ, as it can be found in the reference 
mentioned above. This is an operator that maps kernel- 
like functions, call them w(n, s) where n £ 5' d_1 and 
s £ R, into (filter-like) functions of x £ R d , according to 


Rtw ( x ) = TcJ=TT [ w ( n , ( x |n)) dn, (B2) 
P I JneSd- 1 

Finally, we define lol, the linear transformed kernel of to 
through L, according to the identity 

R^wl((x)) = R^w(L(x)). (B3) 

If uj is the kernel for the filter F, i.e., R^w(x) = E(x), 
then ojl is the kernel for F(Lx). This is precisely the sit¬ 
uation encountered for the squeezed filter provided that 
L is the squeezing transformation. The following exact 
result holds. 

Theorem B.l. (Behaviour under linear transformation) 
Given the filter w(n, s) and an invertible linear transfor¬ 
mation L, then the linear transformed kernel n', s) is 
given by 

wl( n', s) = |detJ4(n)| -1 w(n, |i T n| s). (B4) 

We have here put n' = </>(n) = which is invertible 

on the semisphere, and is its Jacobian matrix. Note 
that the substitution n = <f>~ 1 { n') is intended in the for¬ 
mula above. 

Proof. The application of the adjoint operator to the 
identity that defines col(x) (B3) leads to 

RMix) = 1 [ w(n, (Lx|n))dn. (B5) 

I* I IneS 1 *- 1 

Considering the identity (Lx| n) = <[x| L T n) with L T 
denoting the transposed transformation, it is natural to 
change the variable of the integration to n' = <p(n) = 
, n' £ 5 d_1 , which leads to 


R^wl(x) = R^w(Lx) 

1^ 1 in 'eS^- 1 


w(n, |L T n| (x| n')) 


dn 1 

|detJ4(n)| ’ 
(B6) 


which is the adjoint Radon transform of the function 
u>L(n',s). The statement of Theorem B.l directly fol¬ 
lows by putting (x| n') = s. □ 


For our application of the theorem, we consider a phase 
plane, R 2 , and the linear transformation L corresponds 
to a squeezing transformation (the case of dilation, i.e. L 
is a multiple of the identity is trivial). We write n(9) = 
(cos 9, sin 9) and recover the notation used for the kernel 
in the main text, i.e., co{9, s) = w(n(0), s). The following 
corollary holds. 

Corollary B.l. The general filter that can be obtained 
from a rotational invariant one, denoted by w(s), by di¬ 
lation and rescaling is of the form 

Ua ’ XS) = u 2 (X,9) Wa (uCM)) ’ (B7) 

with uj a (s) defined in Eq. (30) so that ui a ,x{9,s) is nor¬ 
malized as well. X is the rescaling parameter. 

Proof. We parametrize the rescaling transformation by 
We also introduce the the follow¬ 
ing parametrizations n (9) = (cos 9, sin 9) and n'(9') = 
(cos O', sin 9'). Then one can write 

tan#' = A -2 tan 6 
|L T n| = u^ 1 (9 l , A) 

0 ») 

and by substitution the proof is concluded. Note the 
simplification in notation uj a ,x(9, s ) = oj a ,L(n 

At this point, we have established the general results 
on the filter back-projection that we have used through¬ 
out the application of the Inverse Radon method. We 
now turn to the particular problem of the determinating 
the best distribution of cuts C(9), (see Eq. (37)), for the 
Monte Carlo quadrature procedure. 



1. Monte Carlo quadrature. 


By definition, the mean value of the test function and 
A a (e) do not depend on C(9) nor on A [see Eqs. (44) and 
(35)]. Nonetheless, the variance of the estimate (F„ t \) 
has a strong dependence on C(9), that is 


7 Fa. x l C ’ e > M ] ~ 

= f-i JrPa(M— 


x, c (0.s)-(F a , A ))‘ 
M—l 


dsdO 


= Wa^( r+f de \ _ (Fg.N 2 
— M-l\J-% C(8)u 4 (X,8)) M—l ’ 


(B9) 


where the quantity W a ,e = ( f R RIF (t) “*^ dtj has 
been defined. Note that this quantity only relates to 
the unsqueezed distribution and all the A-dependence 
of the variance cr 2 F is contained in the integral 

( f-z c(8)u 4 (\ 8) ') ’ which also depends upon the distri¬ 
bution C. Since our present concern is a suitable choice 
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of C, the integral in the expression above is the main fo¬ 
cus of the analysis. In other words, for fixed a, A, e and 
iff, we wish to determine the distribution C that min¬ 
imises (/J t "| c(e)u*(\, 9 ) ) and a F a , x [ C: > e ’ M } accordingly. 
We can now state the following exact result concerning 
the relevant ratio in Eq. (40). 

Corollary B.2. Assume the state to be described by a 
squeezed distribution W x (x,p), where for X = 1 the distri¬ 
bution is rotationally invariant. Then the optimal distri¬ 
bution of cuts C 0 , which minimises the relevant ratio (40), 
only depends on A through 


C o {0) 


1 

7TU 2 (A, 0) 


(BIO) 


Moreover, for this choice of phase cut distribution the 
value of the relevant ratio does not depend on A and it is 
thus equal to the one for the rotationally invariant state 

(A = 1). 

Proof. We have already pointed out that the dependence 
on the distribution of cuts in the ratio — , is only 

through the integral ^ c(6)u*(\ e) ) appearing in the 

variance. We are thus considering a constrained optimi¬ 
sation problem and can define the Lagrangian functional 


which is positive for any variation SC. C 0 thus gives the 
minimum. Finally, inserting C 0 into <j 2 f [C',e,M] shows 
that the resulting variance Op [C;e, M] is independent 
of A This completes the proof of Colloray B.2. 

□ 


The naive choice of uniformly distributed cuts C = A 
leads to a poor variance for a general squeezed state. 
This is easily seen by straightforward computation of the 
integral 


r+M2 dQ 


r+ir/2 ^4 


da 


I —7r/2 w 4 (M) J —7t/ 2 COS 4 a ( 1 + A 4 tan 2 a ) 5 


+ °° A 2 + A 2 x 2 7T , 2 ,_ 2 \ 

77T 2 \2 dx=-(X- + X 2 ), 

-oo (l + X 2 ) z, 


(B15) 

A" 2 ), 

(B16) 


so that 


_ 2 _ W 0 , e (A 2 + A -2 ) (F aiX ) 2 

<Jf “ x M -1 2 


M - 1 ' 


(B17) 


The first term is dominant for A > 1 (A <C 1), which 
gives the scaling of ~ A 2 (~ A -2 ). 


T£\C,p\ 



C(0)u 4 (\,0) 


p,C(0) dO. 


(BH) 


together with the normalisation condition, C(6)d0 = 
1, that determines the Lagrange multiplier Let 

us consider a small variation of the distribution of cuts, 
C(0) —> C(0)+SC(0). Then, the stationary point equation 
is given by 

0=^m=£; {- 

(B12) 

where ST£ [C, p\ is the linear part of Jz? [C+8C, p] — [C, p] 

with respect to SC. The solution of Eq. (B12) is the sta¬ 
tionary distribution of cuts C a = s ) ■ By integra¬ 

tion, 


r-t-7r/2 


l-n/2 


dO 


r+ir/2 y2 


da 


! (A ,9) 


/_ w /2 cos2 a 1 + A 4 tan 2 a 
r +oc dx 


(B13) 


we find the normalization condition p = 7r 2 . We can 
now check that the stationary distribution C 0 gives the 
minimum of the Lagrangian, by considering the second 
variation around C a , that is 

r +f 

5 2 ^[C 0 ,p] = 7r 3 W a / 5C 2 {0)u 2 {\ 1 9)d0 > 0, 

(B14) 


2. Finite number of cuts. 

Here, we consider the details of the inverse Radon 
method with a finite number of cuts in phase space. As 
remarked in the main text, we assume that the states are 
described by a symmetric W x (x,p) phase space distribu¬ 
tion, where for A = 1 this is the rotationally symmetric 
single photon state. We wish to classify the best way to 
carry out the approximation 

(Fa, a) - W i9a,x(0j) (B18) 

{6j} 

with only a limited number of phase cuts. By the sym¬ 
metry assumed, the mean value of the test function is 
(see Eqs. (31,43,44)) 

= l /_? ^b) (.«<“) 

1 /' + 2 1 

= dO 9a,\=i • (B19) 

Here, it is important that g a , a=i does not depend on 
the phase. The problem thus concerns how to carry out 
the integration of u2 (g \) numerically. In other words, 
if g a X —i does not contain any experimental imprecision, 
the only error associated to the estimate (F a .x) comes 
from the discretisation of f x ^ ■ This problem can 
be addressed with use of the composite middle-point for¬ 
mula. The numerical quadrature is obtained by choos¬ 
ing a partition of the interval [—f, f] into non-uniform 
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m 



FIG. 8 (colour online). Optimal postion for 12 cuts in the 
middle-point quadrature. Here A = 5 and the precision of the 
numerical quadrature is within the 3% of error. 


subintervals, in which the value of the integral is approx¬ 
imated by the product of the length of the subinterval 
times the value of the integrand in its middle point. We 
exploit the parity of the integrand and choose partition 
points symmetric around 9 = 0 , the latter included. The 
complete integral is then approximated as 


1 
7 T 



1 

w 2 (A, 9) 


1 y, AOj 

^« 2 &,A)’ 


(B20) 


where 9j and A 9j are the middle point and the size of the 
j -th subinterval, respectively. The middle points should 
correspond to the cuts where measurements are carried 
out. In analogy with the elementary test, in is the total 
number of cuts over the entire interval [— ^ ]. By com- 

parison with Eq. (B18), the weights are Wj = The 
evaluation at the middle point is a convenient choice be¬ 
cause the integrand is monotonically decreasing (increas¬ 
ing) within the interval [0, + |] ([— f,0]). The middle- 
point evaluation thus provides a trade off between the 
surplus in one half of the interval and the deficiency in 
the other. 

We define the numerical error to be 


g = 



1 y, A 9j 

K fri u 2 (9j,X) 


(B21) 


For fixed A and fixed m, we have numerically optimised 
S with respect to the position of the cuts. As an exam¬ 
ple, we show in Fig. 8 the result for A = 5 and m = 12. 
We observe a concentration of cuts near 9 = 0 , in anal¬ 
ogy with the elementary test. Moreover, as shown in 
Fig. 9, $ decreases roughly exponentially as a function 
of m. The error can, in principle be made arbitrarily 
small by increasing the number of cuts. The bigger A 
is, however, the more cuts are needed to push the error 


FIG. 9 (colour online). The numerical error if as a function 
of the number of cuts m for various values of the squeezing 
parameter A. The plot has been truncated at the onset of 
numerical instabilities that causes the error to oscillate around 
1(T 5 . Notice the nearly exponential decrease for each curve 
and that for larger A more cuts are needed to reach a given 
error. 


below a given threshold. Hence, in an actual experiment 
with a squeezed state, it will be necessary to use a large 
number of cuts to constrain the error from the numeri¬ 
cal quadrature. Also, one would have to justify using a 
highly asymmetric distribution of cuts to estimate an in¬ 
tegral over all angles. This can be argued from the asym¬ 
metry of the filter, which is the main reason for strong 
angular dependence, but such an argument would still 
leave open the question of whether the state itself have 
rapid angular variations in the region with a less dense 
distribution of cuts. Alternatively, one would have to re¬ 
sort to very small intervals between the cuts for all angles. 
If one is using the same number of measurements per cut, 
this would be very costly in terms of measurements since 
a large number of measurements are distributed in a re¬ 
gion with little dependence on the angle. One can to some 
degree compensate for this with an asymmetric distribu¬ 
tion of the number of measurement (see below), but still 
the issue of the precision of the numerical quadrature is 
an uncomfortable question, which it would be desirable 
to avoid, e.g ., by using the elementary test. 

Note also that, for simplicity, the discussion above has 
been carried out by assuming exact knowledge of g a ,\- In 
an actual experiment, where we want to present a proof 
of non-classicality with as few assumptions as possible, 
we do not want to rely on a specific theoretical model for 
the underlying distribution. Instead, one would have to 
justify a specific bound on the error coming from having a 
finite number of cuts by considering the actual measured 
experimental data. This is significantly more complicated 
than the analysis performed here, and we shall not discuss 
it. 

In addition to the error from the finite number of cuts, 
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we also we have to consider the statistical error. If the 
cuts are considered to be statistically independent, the 
variance of the estimate (F a x) is 


a F a 


= £ 


w? 


{6j} 


M, - 1 


(B22) 


where for brevity we have put Wj = satisfying 

Wj = 1. Note that the A dependence is implicit in 
the optimal {9j} and Wj. Also, the variance of g a {9j) is 
independent of 9j. For uniform weights, it is desirable to 
have the same number of measurements Mj for each cut. 
For the optimal quadrature, where Wj are not uniform, 
this is no longer true, since a small weight W 3 allows to 
tolerate a rather big imprecision in g a (9j). Given the to¬ 
tal number of measurement M and m cuts, the optimal 
distribution of measurements per cut Mj can be obtained 
analytically. Formally, the problem consists in the opti¬ 
misation of the quantity in Eq. (B22) with respect to Mj, 
constrained by ^2 3 {Mj — 1) = M — m. The Lagrangian 
for this optimisation problem is 


w: 


■&(Mj,n) = J]~[ +/*E( M i - x )> ( B23 ) 


3=1 


3 = 1 


which has the optimal solution 


Mi- 1 = 


^ Wj 

Vp 


P = 


( ! >: » V 

^ M—m ) 


(B24) 


By eliminating the Lagrange multiplier g, we find (Mj — 
1 ) _ (.m "thU anc [ m i n i ma l value of the variance is 


, - 


(s;.i w, 

M — m 


-** = 


M — m 


(B25) 


The last equality follows from Y^{e } Wj = 1 and estab¬ 
lishes the independence on the squeezing, similar to the 
results obtained for the other methods. Notice that we 
have considered Mj £ M for simplicity; the fact that they 
are natural numbers introduce a small deviation from the 
optimum determined here. 

The total error associated with the estimate (F (h x) is 
given by the sum of all contributions 

<Jp a A + @ + A Q (e). (B26) 

Similar to the other test this can be made independent 
of A provided that m is big enough (see also Fig. 9 and 
discussion; notice although that a large m can decrease 
the strength of the test if M is not very big). Due to 
the large number of ambiguities which would have to be 
resorted in order to make a sensible comparison with the 
elementary test, we will, however, not make a more de¬ 
tailed comparison. 
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